#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <iostream>
#include "WriteMesh.hpp"

void make_mesh(Eigen::MatrixXd &V, Eigen::MatrixXi &T) {
  int x_grid = 2;
  int y_grid = 2;
  int z_grid = 2;
  V.resize((x_grid + 1) * (y_grid + 1) * (z_grid + 1), 3);
  T.resize(x_grid * y_grid * z_grid, 8);
  int index = 0;
  for (int i = 0; i <= z_grid; i++) {
    for (int j = 0; j <= y_grid; j++) {
      for (int k = 0; k <= x_grid; k++) {
        // std::cout<<"i,j,k: " << i << " " << j << " " << k << " " << std::endl;
        V.row(index) =
            Eigen::Vector3d{-1 + 1.0 * k / x_grid, -1 + 1.0 * j / y_grid,
                            -1 + 1.0 * i / z_grid};
        index++;
      }
    }
  }
  index = 0;
  for (int i = 0; i < z_grid; i++) {
    for (int j = 0; j < y_grid; j++) {
      for (int k = 0; k < x_grid; k++) {
        int tmp = k+(x_grid+1)*j+(x_grid+1)*(y_grid+1)*i;
        T(index,0) = tmp;
        T(index,1) = tmp+1;
        T(index,2) = tmp+1+x_grid+1;
        T(index,3) = tmp+x_grid+1;
        tmp = tmp + (x_grid+1)*(y_grid+1);
        T(index,4) = tmp;
        T(index,5) = tmp+1;
        T(index,6) = tmp+1+x_grid+1;
        T(index,7) = tmp+x_grid+1;
        index++;
      }
    }
  }
}

int main(){
    Eigen::MatrixXd V;
    Eigen::MatrixXi T;
    make_mesh(V,T);
    writeHexMeshToVtk("solution_of_sphere.vtk", V, T);
    std::cout << T << std::endl;
}